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ABSTRACT 

For line-driven winds from hot, luminous OB stars, we examine the subtle but important 
role of diffuse, scattered radiation in determining both the topology of steady-state solutions 
and intrinsic variability in the transonic wind base. We use a smooth source function formal¬ 
ism to obtain nonlocal, integral expressions for the direct and diffuse components of the line- 
force that account for deviations from the usual localized, Sobolev forms. As the scattering 
source function is reduced, we find the solution topology in the transonic region transitions 
from X-type, with a unique wind solution, to a nodal type, characterized by a degenerate 
family of solutions. 

Specifically, in the idealized case of an optically thin source function and a uniformly 
bright stellar disk, the unique X-type solution proves to be a stable attractor to which time- 
dependent numerical radiation-hydrodynamical simulations relax. But in models where the 
scattering strength is only modestly reduced, the topology instead turns nodal, with the as¬ 
sociated solution degeneracy now manifest by intrinsic structure and variability that persist 
down to the photospheric wind base. This highlights the potentially crucial role of diffuse 
radiation for the dynamics and variability of line-driven winds, and seriously challenges the 
use of Sobolev theory in the transonic wind region. Since such Sobolev-based models are 
commonly used in broad applications like stellar evolution and feedback, this prompts devel¬ 
opment of new wind models, not only for further quantifying the intrinsic variability found 
here, but also for computing new theoretical predictions of global properties like velocity laws 
and mass-loss rates. 

Key words: radiation: dynamics - hydrodynamics - instabilities - stars: early-type - stars: 
mass-loss - stars: winds, outflows 


1 INTRODUCTION 

For massive, hot stars of spectral types OBA, scattering and ab¬ 
sorption in spectral lines transfer momentum from the star’s intense 
radiation field to the plasma, and so provide the force necessary to 
overcome gravity and drive a stellar wind outflow (see |Puls, Vink| 
& Najarro 2008, for an extensive review). The first quantitative de¬ 
scription of such line-driving was given in the seminal paper by 
jCastor, Abbott & Klein) ( | 1975] ): hereafter “CAK”. Like most wind 
models to date, CAK used the so-called Sobolev approximation 
jSobolevll960f to compute the radiative acceleration. This assumes 
that basic hydrodynamic flow quantitieQ are constant over a few 
Sobolev lengths £sob = v t h/(dv/dr) (for radial velocity gradient 
dv/dr and ion thermal speed vth), allowing then for a local treat¬ 
ment of the line radiative transfer. 

Such a Sobolev approach ignores the strong Tine deshadow- 
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1 Or more specifically, occupation number densities and source functions. 


ing instability’ (LDI) that occurs on scales near and below the 


Sobolev length (Owocki & Rybicki|l984 

Owocki, Castor & Ry- 

bicki 

1988,|Feldmeier, Puls & Pauldrach 

19971 Dessart & Owocki 

2003 

Sundqvist & Owocki 2013 ); this 

is generally believed to 


be the origin for extensive small-scale structure and ‘clumping’ in 


the winds from hot stars (e. g., [Eversberg, Lepine & Moffat|199~8] 
|Puls et al.|2006||Sundqvist, Puls & Feldmeier|2010| . Indeed, while 
Sobolev-based models have had considerable success in explaining 


many qualitative features of stellar winds from OBA stars ( |Paul-| 


drach et al. 1994, 1 

Vink, de Koter & Lamers|2000), recent obser- 

vational analyses (] 

Najarro, Hanson & PulspOll, Sundqvist et al. 

2011 Surlan et al. 

2013[ 

Cohen et al. 2014, Sundqvist, Puls & 

Owocki 2014, Rauw et al. 

2015) that account properly for clump- 


ing systematically infer mass-loss rates in Galactic O-star winds 
that are factors of at least ~ 2 — 3 lower than predicted by such 
Sobolev-based theoretical models. Such mass-loss reductions have 
dramatic consequences for predictions of stellar evolution and feed¬ 
back models (e.g., |Smith|2014) . 

But a further, arguably related shortcoming of the Sobolev ap- 
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proach regards its questionable applicability in regions near and 
below the wind sonic point, where the characteristic scale height 
for variations if =a \p/(dp/dr)\ ~ \v/(dv/dr)\ < -feob- While 
early explorations using full co-moving frame transfer solutions to 
compute the radiative acceleration seemed to suggest quite good 
agreement with the Sobolev approach ( [Pauldrach, Puls & Kudritzki| 
|1986] >, more modern calculations indicate markedly lower radiative 
acceleration in the subsonic region leading up to the wind sonic 
point ( |Lucy||2007a|b[ |2010[ |Krticka & Kubat||2010| |Bouret et al.| 
|2012[ Sundqvist et al., in prep). Such a reduced acceleration in the 
transonic region can have a significant effect on the wind dynam¬ 
ics, perhaps leading to a lower mass-loss rate than predicted by 
corresponding Sobolev models, and thus to better agreement with 
empirically inferred rates. 

As analyzed by Owocki & Puls (1999), the net line-driving in 
the transonic region involves contributions from both a direct com¬ 
ponent of radiative momentum flux absorbed from the underlying 
hydrostatic atmosphere, and a diffuse component of radiation scat¬ 
tered within the local line resonanc^] In the supersonic portions 
of an idealized, steady-state, unclumped and monotonically accel¬ 
erated wind, application of a Sobolev line-transfer treatment shows 
that the escape of radiation scattered within the local Sobolev res¬ 
onance layer attains a fore-aft symmetry, implying then that the 
diffuse component of the line force vanishes in a Sobolev approach 
( |Castor|1974| ). But, as detailed further below, in the transonic re¬ 
gion, fore-aft asymmetries in this escape probability lead to a non¬ 
zero diffuse line-force. 

Nonetheless, because this diffuse term plays little role in driv¬ 
ing the bulk of the supersonic wind outflow, initial dynamical 
explorations of relaxing the Sobolev approach assumed a pure- 
absorption model that completely neglects this diffuse radiation. 
Building on the linear stability analysis of |Owocki & Rybicki] 
< |1984| ), [Owocki, Castor & Rybicki| ( |1988| > used a pure-absorption 
approach in time-dependent simulations to follow the nonlinear 
growth of instabilities on scales near and below the Sobolev length, 
evaluating the direct line-force through non-local, integral calcula¬ 
tions for line-absorption of the underlying stellar continuum. 

|Poe, Owocki & Castor| ( [1990[ hereafter POC) then applied this 
pure-absorption approach to study the nature of associated steady - 
state wind solutions. A major surprise of this analysis was that, 
instead of the unique transonic solution of the usual X-type critical 
point, the solutions near the sonic point of such an absorption-line- 
driven wind exhibit a nodal topology ( |Holzer||1977| ), with both a 
well-defined, steeper-sloped solution, and a degenerate family of 
shallower but still positively sloped solutions (see figure [lj. More¬ 
over, for realistic values of the ratio between the ion thermal speed 
and the sound speed (i.e., v t h/a < 0.3), POC showed that phys¬ 
ically realistic boundary conditions favor the lower-slope, degen¬ 
erate solutions, and that in time-dependent simulations this leads 
to an intrinsic variability that extends down to the subsonic wind 
base. 

Within the context of steady-state models of the transonic 
flow, |Lucy| |2007a|b| ) used non-Sobolev, Monte-Carlo scattering 

2 We emphasize that the “diffuse” component here refers to radiation scat¬ 
tering within the resonance zone for a single, potentially isolated line. This 
is distinct from the multi-line scattered radiation that results when two or 
more distinct lines have overlapping resonances throughout the global wind. 
The latter is commonly modeled with Monte Carlo methods (e.g., |Vink, de| 
Koter & Lamers 2000), using Sobolev-based solutions for location of each 
line resonance, neglecting any net recoil from asymmetries in the radiative 
escape from within that resonance zone. 


line-transfer in a large number (> 10, 000) of spectral lines to com¬ 
pute the total (direct + diffuse) line force. In deriving solutions for 
the associated mass flux through the sonic point, the analytic pa¬ 
rameterization for this numerical, Monte Carlo line force was as¬ 
sumed to depend not on the velocity gradient, as it does in Sobolev- 
based models, but solely on the flow velocity itself, without any 
explicit dependence on the spatial position. As discussed further 
below (see the end of §2.2| ), such an assumed pure-velocity scal¬ 
ing recovers the X-type topology, with one unique positive slope 
solution velocity solution through the sonic point. 

To incorporporate scattering effects in time-dependent dynam¬ 
ical simulations, Owocki (1991) (see also Owocki & Puls 1996 
hereafter OP96) introduced a “smooth source function” (SSF) 
method; this accounts for non-Sobolev asymmetries in the fre¬ 
quency and angle-dependent escape probabilities, but assumes that 
the scattering source function - as a frequency and angle averaged 
quantity - remains smooth and so effectively constant through the 
resonance zone. In this SSF approach, asymmetries in the escape 
near the sonic point lead to a net positive diffuse line-force that 
tends to compensate for the photospheric-line-shadowing reduction 
in the direct component, giving then a total line-force that tends 
to be quite close to the CAK/Sobolev value. Remarkably, for the 
simple case of an optically thin scattering source function with¬ 
out any limb darkening, the base of time-dependent SSF simula¬ 
tions relax to a steady state that matches quite well the standard 
CAK/Sobolev solutions, with now a completely stable transonic 
wind base ( |Owocki|19911|Owocki & Puls|1999| >. 

However, more recent simulations by [Sundqvist & OwocH] 
( [20131 hereafter SO 13) that account for photospheric limb darken¬ 
ing now again show intrinsic base variability and associated struc¬ 
ture in near photospheric layers, in accordance with empirical re¬ 
sults derived from a multitude of O-star observations in various 


wavebands (e.g., Eversberg, Lepine & Moffat 11998, 

Puls et al. 

2006, 

Najarro, Hanson & Puls 2011; Cohen et al. 2011 

Sundqvist 

et al. 

2011, Bouret et al. 2012[ Cohen et al. ;2014). 

This sug- 


gests that even a relatively minor associated reduction in scattering 
strength can have potentially important implications for formation 
of wind structure, and thus also for the empirical derivation of wind 
mass loss rates. This motivates here a renewed examination of the 
dynamics of line driving from the transonic wind base. 

In particular, to clarify and quantify the subtle, but potentially 
crucial role of diffuse, scattered radiation in the initiation and sta¬ 
bility of a line-driven wind through this transonic region, ^extends 
the POC pure-absorption solution topology analysis to account for 
scattering effects within the SSF formalism. A key result is that, 
below some minimum level of scattering, the transonic solution 
topology transitions from the well-defined X-type to the degener¬ 
ate nodal type. Comparison with time-dependent simulations in ^3] 
shows a remarkably close correspondence between the conditions 
for this steady-state topology and the transition from stable relax¬ 
ation to intrinsic base variability. The discussion in ^interprets 
these results in the broader context of modeling line-driven mass 
loss, and outlines some directions for future work. 


2 STEADY-STATE CRITICAL POINT ANALYSIS 

2.1 Direct and Diffuse Components of the SSF Line-Force 

For an isothermal, spherically symmetric, line-driven stellar wind, 
the steady-state equation of motion for variation of radial velocity 
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v with radius r can be written as 

a 2 \ dv 
v J dr 


v — 


GM e fi 2 a 

5 \ Q lines “r 

r z r 


—, (i) 


the line acceleration can be written in terms of the difference be¬ 
tween this term and the inwardly integrated (—) optical depth £_, 
defined by OP96 equation (68), 


where the terms with isothermal sound speed a account for the ef¬ 
fects of gas pressure. Here the gravity GM e ^/r 2 scales with an 
effective stellar mass M e ff = M( 1 — T e ), reduced from the out¬ 
ward force of electron scattering by a fraction set by the Eddington 
parameter T e = n e L / (AttGMc), with the electron scattering 
opacity, M and L the stellar mass and luminosity, G the gravita¬ 
tion constant, and c the speed of light. Since the radiative flux and 
thus line acceleration gimes have a similar inverse-square depen¬ 
dence on radius r, it is convenient to also scale this by the effective 
gravity, 

Flines = GMeff/r 2 ~~ Fdir + Pdiff ’ (2) 

where the latter equality identifies this total line acceleration as 
given by the sum of distinct direct and diffuse components. The 
steady-state topology analysis in POC included only the direct 
component, but we now wish to extend the pure-absorption treat¬ 
ment there to account for scattering terms within the SSF formalism 
dOwocki|199ll , OP96). 

Equations (65) and (67) of OP96 provide relevant expres¬ 
sions for these direct and diffuse forces in terms of integrals over 
thermally scaled frequency displacement from line center x (= 
(y—v 0 )/ Az^d, where Az^d = v 0 v th/c is the Doppler width associ¬ 
ated with ion thermal speed u t h), and over ray index y (= ( p/R ) 2 , 
where p is the ray impact parameter on the projected stellar disk 
of radius R). Following POC, we here compute these integrals in 
terms of associated quadrature sums over discrete grids in x and 
y. For the direct, absorption component of gravitationally scaled 
line-acceleration, we thus write 

Tdir — Tthin ^ ^ ^xy 0 (x 9y^) 5 (5) 

x,y 

with w xy “ w x w y a suitable set of quadrature weights in x and 
y, and with a the power-law index in the assumed line-distribution 
function (see, e.g., equation 52 in OP96). Here <p is the line profile 
function (taken to be a Gaussian), with u = v/v t h the thermally 
scaled radial flow speed, projected along the ray by the direction 
cosine 

(4) 

As in OP96, the line acceleration here is normalized by the opti¬ 
cally thin form for a line of fixed opacity k 0 , 

-p _ ft'o^thT _ ^o^th T e 

thin = 4 tt GMeffC 2 = ^r^TV { 

whereiij^Jthe frequency-integrated line strength K 0 v t h/n e c ~ 1000 
for a typical Galactic O-star |Puls, Springmann & Lennon|2000| . 

The term in equation ([3j represents an associated out¬ 
wardly (+) integrated line optical depth, as given, e.g., by equation 
(66) of OP96. Within the SSF formalism, the diffuse component of 


3 In the notation of Gayley. (19951, the line normalization here is given 
by fto^th /n e c = \QQo a /r(^) ^ where T(a) is the com¬ 

plete Gamma function. The quoted numerical value t hen follows from 
Q « Q 0 « 2000 and a & 0.65 (Gayley|[l995[ |Puls, Springmann &| 
|Lennon|2000) . Note further that in this implementation of the OP96 nota¬ 
tion, the frequency weights w x here must sum to T(a). 


Tdiff = s(r) Tthin ^2 Wx v $ ( x ~ [^- a “ ^+ a ] • ( 6 ) 

x,y 


The t± integral forms (66) and (68) in OP96 are obtained by ap¬ 
plying appropriate boundary conditions to radial integrations of the 
differential form for these outward/inward optical depths, 

dt±(±x,y,r) K 0 p(r)<l>(x - n v (r)u(r)) 

dr ti v {r) 


a relation that will prove useful in the topology analysis below. The 
scattering factor s represents the relative strength of the scattering 
source function to the radial flux. Following OP96, we write this as 


0.5 S 
1 H" fd 1 (r ) *Sthin 


( 8 ) 


with source function S and S = Sthin — (1 — pi(r))/2 for the 
simplest case of an optically thin source function from a star with¬ 
out any limb darkening. 


2.2 Sonic Point Topology in SSF Model 

These SSF forms for the line-acceleration provide the basis for gen¬ 
eralizing the pure-absorption solution topology analysis of POC to 
now account for the effects of the diffuse, scattered radiation. In 
analogy to equation (1) of POC, let us rewrite the equation of mo¬ 
tion 0 in the form 

r 2 dv _ Tlines - 1 + 2a 2 r/GM eff _ N 
GM eff dr ~ V- a?/v ~ D ' U 

As in POC, our analysis here gives the line acceleration in terms 
of both the local radius r and local speed v, allowing one then to 
derive the topology of solutions near the critical, sonic point, where 

N(r c ) D(r c ) = 0 => v c = v(r c ) = a = u c v t h , (10) 


with associated definitions 

_ GMeff 

? 9c — 2 

r 2 

Using L’Hopital’s rule [N/D] c 



; (pc (x) = <p(x - fiycUc) . (11) 

= [dN/dr] c /[dD/dr] c , we obtain 


, 2 9c dN 

Vc 2 dr 


a_ gc dT ii 

o ' rv 


a_ , 9c 

r 2 ^ 2 
A + By' . 


dr 

<9T lines 


dr 


+ v[ 


dr lines 


dv 


(12) 


Applying the above SSF model for the total radiative acceleration, 
along with equation 0 to evaluate the radial derivatives of the opti¬ 
cal depths t±, the coefficients A and B now take the explicit forms, 


































4 Sundqvist & Owocki 



r-r c 


A=-2; B=3 


fm 



W/L 


r-r c 


Figure 1. Illustrations of solution topologies near the critical (sonic) point. 
The left panel shows an X-type case with A = 1, B = 0 and so v' c ± = ±1 
(blue/red lines). The right panel shows a nodal case with A = — 2, B = 3, 
and so v' c+ = +2 (blue line) and v' c _ = +1 (red line). 


A- 


_ 9c <9rii nes /'i 

= 2 c 

= Tthin ^ ] VJ X y { 

x,y 

Kopc ,2 r, — 1 — a /» — 1 — a , 1 — a\l 

- a-0 C [t c+ - s c (t c + + t c _ )J 

Pyc 

~ Sc &c [t~+ ~ t~-] 

+ 2 p yc U c (X - llycU c )(j) c [t~+ - S c (t~+ - t~*j\ } 

(14) 


and 

B 


9c dTi ines 

~2 dv 


^thirijC 

Vth 


^ ] W xy 
x,y 


llyc 0 - Py C U c )(f) c [t c + - Sc(t c + - tc-)] • ( 15 ) 


For pure-absorption (s c — s' c = 0) along a radial ray (y = 0; p y — 
/io = 1; p! y — 0), the last two terms in {l4| l vanish, along with the 
s c elements of the first term, giving then just the POC expression 
for A (cf. their equation 12). Likewise, the vanishing of the s c terms 
in © means that this also recovers the POC expression for B (cf. 
their equation 13). 

From the quadratic equation, we see that solutions for the crit¬ 
ical point slope take the form 


Vc(±) 


B ± V-B 2 + 4 A 
2 


(16) 


As shown by POC, for the pure-absorption case, A < 0 and B > 0. 
If B 2 + 4 A > 0, the critical point is then nodal type, with two 
positive slope solutions, while for B 2 + 4 A < 0), it is focal type, 
with no real solutions for v' c . 

But the addition of the scattering terms now makes it possible, 
for sufficiently large s c , to have A > 0, thus implying a saddle or 
X-type critical point, with two real solutions of v' c that have oppo¬ 
site signs. Figure [I] illustrates the distinct differences between an 
X-type topology (left) vs. a nodal topology (right). 

For nonzero but modest scattering factor s c , the first two terms 
within the sum in 03 are nominally negative, while the final 
term is positive. But the magnitude of the first term is of order 
dgunes/dr « g c /H, where H = a 2 /g c is the atmospheric scale 
height; in comparison, the latter two terms are typically a factor 


H/r c ~ a 2 / g c r c < 0.01 smaller. Moreover, the first term is an 

even bigger factor r 2 /H 2 > 10 4 larger than the gas pressure term 

2/2 
a jr c . 

As such, we may generally neglect all these smaller terms, 
and so approximate the quantity A by just keeping the first term. 
Defining 


y 


Wy_ 

9yc 




X 



/1 + a ’ 

Z c± 


(17) 


we can then write A in the form 

9c 

A ~ ^ Fthin OLKoPc /top &+ 5 (18) 

where we have now defined a ‘topology function’, 

/top = 1 - Sc (i + . (19) 

For the pure-absorption case, with s c = 0, we have / top = 1 and 
so A < 0, giving then a nodal or focal solution topology. To obtain 
A > 0, and thus X-type solutions, f top must become negative, 
which requires some combination of large values for s c and/or for 
the ratio a_/a+. 

Finally, within this context, let us also consider the implied 
topology for the non-Sobolev, Monte-Carlo scattering models of 
|Lucy||2007a|b| >. Since these parametrize the line force as purely a 
function of velocity (thus with nonzero B ~ <9rii nes /c^) but with 
no explicit spatial dependence (and thus vanishing <9rii nes /<9r), 
they have A = a 2 /r 2 c > 0, implying then an X-type solution 
topology, with one unique positive slope solution. Moreover, since 
B ~ g c la ~ a/H a/r c , this positive slope is quite steep, 
«c(+) ~ B ~ a/H , while the negative slope is quite shallow, 
— v c{—) ~ A/B ~ (H/r c )(a/r c ). While this assumed line-force 
parametrization thus gives a unique, non-degenerate solution, fur¬ 
ther work is needed to clarify the physical justification for neglect¬ 
ing any explicit spatial dependence in parametrizing the Monte 
Carlo line force. 


2.3 Sobolev optical depth for forward/backward streams, r± 


Appendix A details a formal second-order Sobolev analysis for es¬ 
timating the expected values of a± and thus A. To build insight, let 
us consider here a simpler, more heuristic analysis for this. Within 
the Sobolev approximation, the optical depth integrals scale in pro¬ 
portion to a locally defined Sobolev optical depth, t ~ r = n 0 p£, 
where t = v t h/v' is the Sobolev length. But for the critical-point 
forward/backward streaming optical depths t c ±, the relevant eval¬ 
uations center on regions that extend roughly a Sobolev length be¬ 
low/above the sonic point; in terms of sonic point values of the 
Sobolev length £ c = v t h/v' c and optical depth r c = n 0 p c ^c, the as¬ 
sociated forward/backward streaming Sobolev optical depths scale 
as 


Z°± = 

T c 2 T c 


lT l“ 

2 Pc 


i±^ = i 

2 a 2 a 


( 20 ) 


where r' c represents the local radial gradient of r c . Here the second 
approximation assumes that the velocity gradient is itself nearly 
constant (i.e., that v'/ « 0), while the latter forms use the constancy 
of mass flux pv for a steady, nearly planar outflow. 

Since t c ± ~ r c ±, we then see that 


a~ 

a + 



l + (l + a) 


Vth 

? 

a 


( 21 ) 


where the latter approximation applies first order expansion in the 
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small parameter v t h/a < 1. Applying this in {19} , we can estimate 
a minimum value for the scattering constant to obtain / top < 0 and 
thus A > 0 for an X-type critical solution, 


Sc > Sx 



1 

2 _|_ (! + o0^th 


( 22 ) 


In the strict Sobolev limit v t h —> 0, obtaining an X-type 
neutral point would require a sonic region scattering factor s c > 
sx = 1/2. For a simple optically thin scattering source-function, 
a surface intensity that is either constant or limb darkened has 
s(r) ^ 1/2 for all r ^ R (see equation [8|. Such a strict Sobolev 
limit model would thus have to have photospheric limb brightening 
to achieve the A > 0 needed for an X-type solution topology. 

But for a small but nonzero thermal speed ratio v t h/a < 1, 
equation ( [22] ) gives sx < 1/2; for example, for a m 0.65 and 
vth/a ~ 0.28, we find sx ~ 0.4, making it easier to get A > 0 
and so a stable, X-type critical solution. This can thus help explain 
why the base of SSF models without limb darkening relax to a sta¬ 
ble steady state ( |Owocki||1991[ |Owocki & Puls||1999| ). However, 
when limb darkening is included, the lower value of s c can again 
lead to A < 0, and so a nodal topology, with its associated intrin¬ 
sic variation at the wind base, as found in the SSF simulations by 
S013. 

Within the more formal second-order Sobolev analysis of Ap¬ 
pendix A, figure | A1 1 illustrates the asymmetries in t c ± that lead to 
a- > a+, while figure |A2| plots the resulting variations of a_/a + 
and sx with Vth/a. 


3 COMPARISON WITH TIME-DEPENDENT SSF 
SIMULATIONS 

Let us now compare and interpret the above analytic results in the 
context of full numerical time-dependent hydrodynamic wind sim¬ 
ulations based on the SSF model for line-driving. 


3.1 Basic Model Description 

The simulations here use the numerical hydrodynamics cod^] 
VH-1 to evolve the conservation equations of mass and momentum 
for a spherically symmetric and isothermal outflow. As in previ¬ 
ous SSF simulations (e.g., |Owocki & Puls|1999[ SO 13), we imple¬ 
ment radiation line-driving following the method described in 
All results presented here adopt the same stellar and wind param¬ 
eters as in SO 13, given by their Table 1, which are typical for an 
O-star in the Galaxy. For the models in j |3.2.1| and §3.2.3| we fix 
the ratio between ion thermal speed and isothermal sound speed to 
vth/a = 0.28, representative of driving by carbon, oxygen, and 
nitrogen atoms (as in SO 13; see also discussion in POC); in §3.2.2| 
we explore different fixed values of v t h/a to study its influence on 
the stability of the wind base. 

Each simulation evolves from a smooth, CAK-like initial con¬ 
dition, computed by relaxing to a steady state a time-dependent 
simulation that uses a CAK/Sobolev form for the line-force. To 
prevent artificial structure due to numerical truncation errors we 
use a fixed evolution time-step of 2.5 s, rather than setting this step 
to a fixed fraction of the courant time (see discussion in POC). As 

4 The VH-1 hydrodynamics computer-code package has been developed 
by J. Blondin and collaborators, and is available for download at: http: 
[7/wonka . physics . ncsu . edu/pub/VH-1/ 


in previous work, the lower boundary at the assumed stellar sur¬ 
face fixes the density to a value ~ 5 — 10 times that at the sonic 
point. Since we are interested here in topology and intrinsic wind 
variability, we do not introduce any explicit base perturbations (see 
SO 13, for examples of SSF simulations including such explicit per¬ 
turbations); this means all structure seen in the models here is truly 
intrinsic to the line-driven wind. The spatial grid uses 500 discrete 
radial mesh-points out to r — 1.5AL, with a step-size that increases 
linearly with radius, to increase the resolution in the transonic re¬ 
gion critical for the analysis here. Finally, all simulations use a sin¬ 
gle ray coordinate at y — 0.5 to approximate the effects on the 
line-force from the finite stellar disc; tests have shown that such 
a single ray typically gives a radiative acceleration that deviates 
< 10% from that computed from full angle integration (see SO 13; 
their Figure A.l). 


3.2 Simulation Results 

3.2.1 Fixed S /Sthin factors 

Figure [2] shows results from SSF simulations, using three simple 
forms S'/S'thin = 1, 0.95, and 0.9 for the source function that sets 
the scattering factor s(r). (See equation]!]) The first of these three 
cases (shown in the upper row of the figure) thus assumes an opti¬ 
cally thin scattering source function from a uniformly bright stel¬ 
lar disk, while the second (middle row) and third (lower row) cases 
simply reduce this source function by respectively 5% and 10%. As 
we now demonstrate, such simple variation of the scattering term 
s allows us to examine carefully the transition between nodal and 
X-type topology in our numerical simulations. 

In the associated simulation results shown in Figure[2] the left¬ 
most panels display the ratio a_ /a+ (equation[l7]), computed every 
1000 s at the mesh-point closest above the critical sonic point r c . 
For a given scattering term s(r c ), the next panels plot the topology 
function / top at this mesh-point. The remaining two panels display 
radius vs. time color-plots of mass-loss rate M = Airr 2 pv and 
density. 

As evident from the figure, after an initial adjustment to the 
new force conditions, the standard simulation with SySthin = 1 
relaxes to a stable wind base, wherein the line-drag of the diffuse 
radiation ( |Lucy|198'4]|Owocki & Rybicki|1985| exactly cancels the 
LDI at the stellar surface; in this case, self-excited structure devel¬ 
ops only near the outer boundary at r « l.bR here. (See SO 13 for 
structured simulations extending to larger radii.) 

Building on the analysis in ^2] we can now interpret the rela¬ 
tive stability of the base in terms of the / top function at the sonic 
point. Even though s(r c ) < 0.5 (reflecting the fact that the sonic 
point is located slightly above the stellar surface), the a_/a+ ratio 
yields / top < 0 at all times for this standard case (see first and 
second panels in top row), allowing then the wind around the sonic 
point to relax to a unique X-type solution. 

But in the runs with S/S t hin = 0.95 and 0.9, the lower values 
of s(r c ) never allow the base to completely stabilize. At certain 
time-s-eps-converted-to.pdf / top again becomes positive, implying 
a nodal wind topology. In the time-dependent simulations here, this 
gives rise to the same type of intrinsic global variability as found 
in the POC pure-absorption models, and in the limb-darkened SSF 
simulations by SO 13. 

Note that the run with S'/S'thin = 0.95 exhibits quite long pe¬ 
riods where the wind appears almost stable, but as / top approaches 
and exceeds unity, the wind reacts by admitting the nodal topol¬ 
ogy and its non-unique shallow-slope solutions. As compared to 
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Figure 2. Time-dependent SSF simulations computed with scattering source functions S/ 5th = 1 (upper row), 0.95 (middle row), and 0.9 (lower row). All 
ordinates display time in ksec. From left to right the panels in each row show the ratio a_/a+, the topology function at the wind sonic point /to P (r c ), and 
color plots of mass-loss rate (log dM/dt = log M [Mq /year]) and density (log p [ g/cm 3 ]), where the abscissae show the radial coordinate in units of the 
stellar radius, r/R*. See text. 


the transitional simulation with S'/S'thin = 0.95, the run with 
S/Sthin = 0.9 has an even lower scattering term. After some ini¬ 
tial adjustment-time, this now results in a strong and almost regular 
variability reaching all the way down to the wind base. Again times 
with /top ^ 1 correspond remarkably well to periods of strong 


variability, in between which the wind initially appears to relax, 
but then fails to settle down to a steady-state. 
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Figure 3. Sonic point topology function for different fixed ratios of ion 
thermal to sound speed. The red solid diamonds show results from SSF 
simulations with S/Sthin = 1 and the black solid line compares this to 
the heuristic analytic prediction eqn. [fusing s(r c ) = 0.47 (appropriate 
for simulations with S/S t hin = L see text). The dashed line marks the 
/top = 0 line where the wind transitions from X-type to nodal topology. 


3.2.2 Dependence on ion thermal speed 

As mentioned, the standard SSF models above all assume v t h/a = 
0.28. But the heuristic analysis in §2.3[ predicts that as the Sobolev 
limit v t h —^ 0 is approached, obtaining a stable X-type solution 
requires a sonic point critical scattering term s c > s x = 1/2. 
And since an optically thin source function has s x < 1/2 at the 
sonic point, this implies that as the ratio v t h /cl is lowered, even 
SSF models with S/S t hm = 1 should eventually become nodal 
and so exhibit an unstable base. 

To investigate this, we ran a set of standard SSF simula¬ 
tions with S/Sthin — 1 but now with different vth/a ratios. Fig¬ 
ure [3] compares the asymptotic, steady-state values of / top for 
these models to the simple analytic scaling equation © (using 
s(r c ) = 0.468, appropriate for our standard optically thin source 
function). The figure shows that while the non-linear SSF simu¬ 
lations have larger escape asymmetries (i.e. more negative values 
of /top) at vth/a = 0.28, the numerical models agree quite well 
with the simple analytic scaling at lower thermal speeds. As an¬ 
ticipated, the simulations approach unstable conditions when the 
thermal speed is lowered, crossing the / top = 0 point that marks 
the transition from X-type to nodal topology at vth/a ~ 0.07. 
Inspections of simulations with v t h/a < 0.07 further show that 
such models indeed exhibit a very chaotic wind base with exten¬ 
sive structure reaching all the way down to the lower boundary, as 
demonstrated by Fig. [4] This figure plots density colour-maps of 
two models with vth/a — 0.07 and 0.05, illustrating how the first 
of these simulations (which has / top very close to zero, see Fig. [3]) 
is marginally stable, whereas the latter displays extensive structure 
already near the wind base. 

3.2.3 Limb-darkened models 

The SSF simulations above clearly demonstrate how scattering in 
the transonic wind region controls the admitted topology through 
the critical wind sonic point region. Building on these results, let us 
thus now examine the somewhat more physically realistic case of 


Figure 4. Density color-maps of two SSF simulations with S/Sthin = 1 
and vth/a = 0-05 (left) and 0.07 (right). See text. 


an Eddington-limb-darkened source function (like that adopted by 
SO 13), returning again to the standard thermal speed with vth/a — 
0.28. 

Figure[5]shows the same set of plots as in Figure[2] but now for 
such a model. As discussed in §2, a limb-darkened source function 
naturally lowers the scattering term s(r c ) and can thus easily reach 
/top ^ 1 and so nodal topology. And, indeed, Figure [5] makes it 
immedietely clear that the strong variability extending down to the 
wind base in such limb-darkened SSF simulations is a direct con¬ 
sequence of /top reaching values above unity during semi-regular 
intervals. As in the simulations above, this then gives rise to nodal 
topology solutions lying on the degenerate shallow-slope branch. 
Moreover, as discussed in detail in SO 13, in such limb-darkened 
SSF simulations the lower scattering term also means the diffuse 
line-drag no longer cancels the LDI at the wind base, leading to 
non-linear feedback and presumably further enhancement of struc¬ 
ture in these near-photospheric layers. 


4 DISCUSSION, CONCLUSIONS, AND OUTLOOK 

The analysis here highlights the critical role played by the diffuse, 
scattered radiation field in the dynamics and variability of line- 
driven winds. 

In a non-Sobolev model, in which the line force depends ex¬ 
plicitly on velocity and radial coordinate (rather than on the veloc¬ 
ity gradient), the sonic point constitutes a “critical point” of vital 
importance for the wind dynamics. By extending here the pure- 
absorption analysis of POC, we show that standard, stable X-type 
wind solutions are admitted only in cases where the combination 
of the scattering term and the fore-aft asymmetry in escape is large 
enough. Within the SSF formalism, this is controlled by the sign 
of a sonic-point “topology function” / top ; a negative / top gives 
X-type topology, whereas a positive / top leads instead to nodal 
topology, with quite different behavior and an allowed branch of 
non-unique wind solutions. 

For the standard case of v t h/a — 0.28 and an idealized opti¬ 
cally thin scattering source function from a uniformly bright stel¬ 
lar disk, SSF simulations give numerically / top < 0, with then 
an associated relaxation to a stable wind base. But in cases where 
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Figure 5. Like Figure[2] but for a limb-darkened source function. 


the source function is only mildly reduced from this idealized case 
(as in simple limb-darkened models), / top becomes positive at cer¬ 
tain times, implying a degenerate nodal topology. In the spheri¬ 
cally symmetric simulations here, this leads to intrinsic structure 
and a semi-regular variability extending down to the wind base. 
This may be crucial for setting the near-photo spheric scale of var¬ 
ious phenomena observed in hot star winds (e.g., wind clumping 
and discrete-absorption components ; |Eversberg, Lepine & M offat 
|1998| |Kaper et aT7|| 1 999| |Puls et aL||2006| ), and follow-up studies 

will examine how this variability manifests and evolves in multi¬ 
dimensional simulations, where the spherical shell structure is bro¬ 
ken up by thin-shell and other hydrodynamic instabilities ([Dessart 
|& 0wocki|2003l|2005] >. 

The analysis here further illustrates quite vividly how very 
sensitive the wind dynamics is to details regarding the scattered ra¬ 
diation field in the transonic region. As such, it seriously questions 
the validity of using Sobolev theory to compute the line force in this 
region, thus challenging the Sobolev-based wind models often used 
in broad applications like stellar evolution (e.g., |Vink, de Koter &| 
Lamers 2000). This motivates development of a new generation of 
models, able to quantitatively compute global wind properties ac¬ 
counting properly also for the diffuse component of the line force. 

Indeed, full co-moving frame radiative transfer calculations 
typically find a significantly reduced source function in the tran¬ 
sonic region ( [Krticka & Kubat||2010[ |Bouret et al.||2012[ Sund- 
vist et al., in prep.), and initial experimentation including a sim¬ 
plified escape-integral source function method (EISF, |Owocki &| 
|Puls|1999} in time-dependent, dynamical simulations shows signif¬ 
icantly lower mass-loss rates and steeper velocity fields than pre¬ 
dicted by Sobolev-based models ( [Owocki & Puls||1999| ). In view 
of the analysis here, it is interesting to note that while such EISF 
models certainly seem to lie on the nodal solution branch, they ap¬ 
pear to relax to the unique steep solution, rather than admitting the 
degenerate shallow-slope solutions found in the SSF simulations 
here. 

Also the time-independent Monte-Carlo scattering models of 
the transonic flow by |Lucy| ( |2010| ), find significantly lower mass- 
loss rates for main-sequence O-stars than predicted by the Monte- 
Carlo Sobolev simulations by |Vink, de Koter & Lamers] j2000| . 
However, in addition to never driving the wind beyond a few times 


the sound speed, the assumed pure-velocity scaling of the line-force 
in the Lucy (2007a b 2010 ) models effectively suppresses appear¬ 
ance of a nodal topology, leaving instead only a single X-type so¬ 
lution (see discussion in §2.2| ). 

Further work, e.g. using co-moving frame transfer to compute 
the radiation force, is thus needed to clarify which solution is ob¬ 
tained in dynamical simulations of line-driven winds, and under 
what circumstances the base of such winds relax to a steady state. 
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Figure Al. Terms from the integrand in the a± integral equation 
plotted vs. the negative of the critical point comoving frequency —x = 
u c — x, for cases with u t h/a = 0.07 (left) and v t h/a = 0.28 (right). 
The blue and red dashed curves show respectively forward/backward factors 
(r c /t c ±) 1+a ; the dotted black curve shows 27r</> 2 (x — u c ); and the filled- 
area curves show the corresponding products of these, representing the full 
integrand of 03 - Note that the red shaded area, representing a_, is larger 
than the blue shaded area, representing a+, especially in the right panel 
with larger u t h /a. The headers in each panel label both the assumed value 
°f ^th / o and the computed value of a_ /a+. 


APPENDIX A: SECOND-ORDER SOBOLEV ANALYSIS 

To complement the simplified forward/backward optical depth es¬ 
timates in §2.3[ let us now carry out a more formal “second-order 
Sobolev” analysis of the full optical depth integrals. For simplicity, 
we assume a purely radial ray with y — 0 and thus p y = 1. From 
equation (66) of OP96, we can write the forward-streaming optical 
depth at the critical point as 


tc+(x) = t+(x, 0, r c ) 


+ At c+ (x), (Al) 

^max 


where the first two terms account for the line-opacity cutoff (taken 
here to have a value ftmax ~ 10 6 ft e ) and the Shuster-Schwarzschild 
lower boundary condition at the stellar radius R; the final term is 
given by the outward radial integral 


1 + — (r — r c ) 
Pc 


1 + 


(j)(x — u(r))dr 
dx 


rr c 

At c +(x) = / fto p{r)4>(x - u(r )) dr 

JR 

rr 

KoPc I 
Jr 

rx 

KoPc l 

Jx-u(R) L a U c 

T c f [2 + — (x — £C)1 Mx) dx 
Jx-u c La J 

Tc [^1 - ( erf ( X ) “ erf ( X - u c )) 

+ ^ ~ UC ^ _ 


(A2) 

(A3) 

(A4) 

(A5) 

(A6) 


Here the intermediate forms assume a nearly linear variation in den¬ 
sity and velocity within an approximately planar, steady-state flow 
near and below the critical (sonic) point. 

The corresponding backward-streaming optical depth takes 
the form (cf. equation 68 of OP96), 


t c -(x) = t- (-x, 0, r c ) = + At c - (x ), (A7) 

^max 
















10 Sundqvist & Owocki 



Figure A2. Left: The ratio a_/a+ plotted vs. v^/a (solid curves), com¬ 
pared with the analytic expansion form a_/a+ ~ 1 + (1 + a)v t h/a 
(dotted; cf. equation |2l) and the simple unit-slope linear form a_/a+ « 
1 + v t h/a (dashed). Right: Associated variations of the minimum source 
function for an X-type topology, as defined by equation (A12) . In both pan¬ 
els, the filled circles highlight results for the left and right panels in figure 


[ATI 


where 

A t c -(x) 


px — u, 

TC J- oo 


To [fl- 


r oo 

(x) = / k 0 p(r)4>(x — u(r)) dr 
Jr c 

2 -f ^2 (x — x)j 4>(x) dx 


(A8) 

(A9) 


These analytic forms for t c ± can be used to compute the terms 
that determine the sign of A, which here depend on the full fre¬ 
quency integral forms of equation GZf 


a± = 


f 


fix) 


dx . 


(All) 


To produce X-type critical solutions with A > 0 , the critical-point 
scattering coefficient s c must now exceed, 

SX = . (A 12) 

+ a + 

Figure |A1| plots the frequency variation of the relevant inte¬ 
grand terms in the integral equation fTT] ) for a±, comparing results 
for vth/a = 0.07 (left) vs. v t h/a = 0.28 (right). Figure [A2| plots 
the ratio a_ /a + plotted vs. vth/a (solid curves), with the filled cir¬ 
cles highlighting results for the two cases shown in figure [AT] the 
dotted and dashed curves compare simple linear forms with slopes 
of unity and 1 + a = 1.65, respectively. 

For the standard case vth/a = 0.28 and a = 0.65, evaluation 
gives a_/a+ ~ 1.35 and sx = 0.425, results that are actually in 
quite good agreement with the more heuristic analysis in §2.3| 






